{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.2 CoefficientFunctions\n",
    "====\n",
    "\n",
    "In NGSolve, `CoefficientFunctions` are representations of functions defined on the computational domain $\\Omega$. Examples are expressions of coordinate variables $x, y, z$ and functions that are  constant on subdomains. Much of the magic behind the seamless integration of NGSolve with python lies in `CoefficientFunctions`. This tutorial introduces you to them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "mesh = Mesh (unit_square.GenerateMesh(maxh=0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ngsolve.fem.CoefficientFunction at 0x7f6e44cfad00>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myfunc = x*(1-x)\n",
    "myfunc   # You have just created a CoefficientFunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ngsolve.fem.CoefficientFunction at 0x7f6e44cfa518>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x        # This is a built-in CoefficientFunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Draw the function\n",
    "\n",
    "Use the `mesh` to visualize the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(myfunc, mesh, \"firstfun\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Evaluate the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.16000000000000003"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mip = mesh(0.2, 0.2)\n",
    "myfunc(mip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `myfunc(0.2,0.3)` does not work: You need to give points in the form of `MappedIntegrationPoint`s like `mip` above. The `mesh` knows how to produce them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examining functions on sets of points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "point=(0.00,0.20), value=0.00000\n",
      "point=(0.10,0.20), value=0.09000\n",
      "point=(0.20,0.20), value=0.16000\n",
      "point=(0.30,0.20), value=0.21000\n",
      "point=(0.40,0.20), value=0.24000\n",
      "point=(0.50,0.20), value=0.25000\n",
      "point=(0.60,0.20), value=0.24000\n",
      "point=(0.70,0.20), value=0.21000\n",
      "point=(0.80,0.20), value=0.16000\n",
      "point=(0.90,0.20), value=0.09000\n"
     ]
    }
   ],
   "source": [
    "pts = [(0.1*i, 0.2) for i in range(10)]\n",
    "vals = [myfunc(mesh(*p)) for p in pts] \n",
    "for p,v in zip(pts, vals):\n",
    "    print(\"point=(%3.2f,%3.2f), value=%6.5f\"\n",
    "         %(p[0], p[1], v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may plot the restriction of the `CoefficientFunction` on a line using matplotlib. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "px = [0.01*i for i in range(100)]\n",
    "vals = [myfunc(mesh(p,0.5)) for p in px]\n",
    "plt.plot(px,vals)\n",
    "plt.xlabel('x')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interpolate a `CoefficientFunction`\n",
    "\n",
    "We may `Set` a `GridFunction` using a `CoefficientFunction`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=1)\n",
    "u = GridFunction(fes)\n",
    "u.Set(myfunc)\n",
    "Draw(u)         # Cf.: Draw(myfunc, mesh, \"firstfun\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The `Set` method interpolates `myfunc` to obtain the grid function `u`.\n",
    "\n",
    "* `Set` does an Oswald-type interpolation as follows:\n",
    "    - It first zeros the grid function;\n",
    "    - It then projects `myfunc` in $L^2$ on each mesh element;\n",
    "    - It then averages dofs on element interfaces for conformity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Integrate a `CoefficientFunction`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can numerically integrate the function using the mesh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.166666666666666"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Integrate(myfunc, mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Differentiate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no facility to directly differentiate a `CoefficientFunction`. But you can interpolate it into a `GridFunction` and then differentiate the `GridFunction`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ngsolve.fem.CoefficientFunction at 0x7f6deaab8518>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u.Set(myfunc)\n",
    "gradu = grad(u)\n",
    "gradu[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(gradu[0], mesh, 'dx_firstfun')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously the accuracy of this process can be improved for smooth functions by using higher order finite element spaces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Vector-valued `CoefficientFunctions`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, `gradu` provided an example of a vector-valued coefficient function. To visualize it, click on `Visual` menu in GUI and check `Draw Surface Vectors`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(gradu, mesh, \"grad_firstfun\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also define vector coefficient expressions directly: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "vecfun = CoefficientFunction((-y, sin(x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw(vecfun, mesh, \"vecfun\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Expression tree\n",
    "\n",
    "Internally, coefficient functions are implemented as an expression tree made from building blocks like `x`, `y`, `sin`, etc., and arithmetic operations.\n",
    "\n",
    "E.g., the expression tree for `myfunc = x*(1-x)` looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coef binary operation '*', real\n",
      "  coef coordinate x, real\n",
      "  coef binary operation '-', real\n",
      "    coef N5ngfem27ConstantCoefficientFunctionE, real\n",
      "    coef coordinate x, real\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(myfunc) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
